function [parms_fitted, gesol, FitVal, FLAG_SUCCESSFUL_FSOLVE] = ...
    FitModel20220926_comp_EU_stdU_stdmu3(parms, b_guess, kappa_guess, log_sigma_x_guess,...
        U_mean_target, U_std_target, ave_mu3_vol_target, NumTries, fsolve_opts)
% target moments:
% E[U], vol[u]
% parameters:
% b, kappa

    if b_guess<=0; error('check b_guess'); end
    if kappa_guess<=0; error('check kappa_guess'); end
    
    fixedpt_tol = 1e-6; % speed things up
    
    N_Nxz = repmat(parms.grid_N(:),[1 parms.Nx parms.Nz]);
    U_Nxz = 1 - N_Nxz;
    x_Nxz = permute(repmat(parms.grid_x(:),[1 parms.NN parms.Nz]),[2 1 3]);
    phi_of_x = @(x) 1./(1 + exp(-x));
    phi_Nxz = phi_of_x(x_Nxz);
    log_phi_Nxz = log(phi_Nxz);
    input_mu3_Nxz = (1 - U_Nxz).*U_Nxz.*(1 - 2*U_Nxz).*(log_phi_Nxz.^3);
    
    FLAG_SUCCESSFUL_FSOLVE = false;
    
    for iter = 1:NumTries
        fprintf('iter %g with b_guess=%g, kappa_guess=%g.\n', ...
            iter, b_guess, kappa_guess);
        [XFit, FitVal, EXITFLAG] = fsolve(@FitFun,[log(b_guess); log(kappa_guess); log_sigma_x_guess],fsolve_opts);
        if EXITFLAG>0
            fprintf('Found roots on iteration %g.\n', iter);
            FLAG_SUCCESSFUL_FSOLVE = true;
            break;
        else
            fprintf('iter %g, fitted vals (not solution):\n b=%g, kappa=%g.\n', ...
                iter, exp(XFit(1)), exp(XFit(2)));
            b_guess = 0.3 + 0.6*rand;
            kappa_guess = 2*rand;
            log_sigma_x_guess = log(0.01 + 0.2*rand);
        end
    end
    
    % do a last compute
    parms_fitted = parms;
    parms_fitted.b(:) = exp(XFit(1));    
    parms_fitted.kappa(:) = exp(XFit(2));
    parms_fitted.sigma_x(:) = exp(XFit(3));

    fprintf('kappa=%g, b=%g, kappa=%g.\n', ...
        parms_fitted.b(1), parms_fitted.kappa(1));

    gesol = SolveFixedPoint_Model_20220926(parms_fitted,fixedpt_tol);
    parms_fitted.FitVal = FitVal;

    if ~FLAG_SUCCESSFUL_FSOLVE
        warning('Model fit: not all moments are satisfied.');
    end

    function Z = FitFun(X)
        
        Z = NaN*zeros(3,1);
        
        parms_temp = parms;
        parms_temp.b(:) = exp(X(1));
        parms_temp.kappa(:) = exp(X(2));
        parms_temp.sigma_x(:) = exp(X(3));

        gesol_temp = SolveFixedPoint_Model_20220926(parms_temp,fixedpt_tol);
        if gesol_temp.VFIerr>fixedpt_tol
            warning('FitModel20220926_ver1: VFI error = %g is greater than tolerance of %g.', ...
                gesol_temp.VFIerr, fixedpt_tol)
            bComputeMoments = false;
        else
            bComputeMoments = true;
        end
        
        %% Newer version
        
        if bComputeMoments

        [Phi_Nxz,ss_dist_Nxz] = MakeSparseTransMatrix_Nxz(parms_temp,gesol_temp);

        % NEW CODE: quarterly average
        [mean_U, std_U] = Calc_quarterly_ave_NHz(Phi_Nxz,ss_dist_Nxz,U_Nxz);

        % cross-sectional mu3 vol
        std_mu3_cons_1y = CalcTimeSeriesDiffVol(input_mu3_Nxz,ss_dist_Nxz,Phi_Nxz,12);
        
        Z(1) = (mean_U - U_mean_target)/U_mean_target;
        Z(2) = (std_U - U_std_target)/U_std_target;
        Z(3) = (std_mu3_cons_1y - ave_mu3_vol_target)/ave_mu3_vol_target;
        
        end
        
    end

end